########################################################################
# PLS 206 — Week 5: Model selection
# Grey Monroe, Fall 2026
#
# Every block is tagged:
#   [CORE]   the analysis: code you would run on your own data
#   [HANDY]  useful, not essential: wrangling, plotting, small helpers
#   [DEMO]   for the lecture only: simulated data, sums by hand, and the
#            code behind slide figures
########################################################################

# [CORE]
# install.packages(c("ggplot2", "MuMIn", "car", "MASS"))
library(ggplot2)
library(MuMIn)   # AICc(), model.sel(), dredge(), sw(), model.avg()
library(car)     # vif()
library(MASS)    # stepAIC()

# ---- Which tool for which question -----------------------------------

# | You want… | Use | |---|---| | A score for one model | AIC(m),
# AICc(m) (MuMIn), BIC(m) | | To compare a few models you wrote down
# before looking | model.sel(m1, m2, m3) | | A test for one nested pair
# | anova(m_small, m_big) | | To find collinear predictors | vif(m)
# (car): drop the worst, refit, repeat | | To search every subset
# (exploration) | options(na.action = "na.fail"), then dredge(full,
# fixed = ~ control) | | How much support each variable has | sw(dd) |
# | Estimates across the well-supported models |
# model.avg(get.models(dd, subset = delta < 2)) | | To judge
# predictions | fit on part of the data, predict the rest |

########################################################################
# PART A. Too many models
########################################################################

# ---- 1) The data: the cover-crop farms -------------------------------

# The cover-crop data from the Week 2 live coding: spring biomass
# (biomass, kg/ha) of the same cover crop on 50 farms, with each farm's
# location, winter weather and soil. clay goes, because sand + silt +
# clay = 100 (Week 2). Simulated: only growing degree days (gdd) drives
# biomass, so we can check what each method finds.

# [CORE]
farms <- read.csv("https://greymonroe.github.io/PLS_206/data/cover_crop.csv")
farms <- subset(farms, select = -c(farm, clay))
head(farms)
nrow(farms)

# [HANDY]
P <- setdiff(names(farms), "biomass")     # the 15 candidate predictors
long <- do.call(rbind, lapply(P, function(v) data.frame(var = v, x = farms[[v]], biomass = farms$biomass)))
long$var <- factor(long$var, levels = P)
ggplot(long, aes(x, biomass / 1000)) +
  geom_point(colour = "grey50", size = 1) +
  facet_wrap(~ var, scales = "free_x", nrow = 3) +
  labs(x = NULL, y = "Biomass (t/ha)") +
  theme_classic()

# ---- 2) How many models could we fit? --------------------------------

# Each predictor is in or out, so 15 predictors give 2^15 models,
# before any interactions or curves.

# [CORE]
ncol(farms) - 1         # candidate predictors
2^15                    # possible models (each one in or out)

# ---- 3) Add predictors: R² only goes up ------------------------------

# Fit on 25 farms, add predictors one at a time, best first (the
# strongest correlation with biomass on those 25 farms), until all 15
# are in. Keep the R² on the 25 farms and the prediction error (RMSE)
# on both halves.

# [DEMO: fit on 25 farms, add predictors best first]
set.seed(5)
tr <- sample(50, 25)                                          # the 25 farms we fit to
ord <- names(sort(-abs(cor(farms[tr, ])[P, "biomass"])))      # predictors, best first
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))
of <- do.call(rbind, lapply(0:15, function(k) {
  fm <- if (k == 0) biomass ~ 1 else reformulate(ord[1:k], "biomass")
  m <- lm(fm, data = farms[tr, ])
  data.frame(k = k, r2 = summary(m)$r.squared,
             train = rmse(farms$biomass[tr], fitted(m)),
             test  = rmse(farms$biomass[-tr], predict(m, farms[-tr, ])))
}))
round(of, 3)    # R-squared 0.57 with 1 predictor, 0.88 with 15

ggplot(of, aes(k, r2)) +
  geom_line(colour = "blue", linewidth = 1) + geom_point(colour = "blue", size = 2) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "Number of predictors in the model", y = "R² on the 25 farms it was fitted to") +
  theme_classic()

# ---- 4) ...but the other 25 farms get worse --------------------------

# The error on the 25 farms the model was fitted to falls from 1.04 to
# 0.56 t/ha. On the other 25 farms it is lowest with one predictor
# (gdd, 1.16 t/ha) and reaches 2.57 t/ha with all 15.

# [DEMO: prediction error on both halves, and predicted vs observed]
ol <- rbind(data.frame(k = of$k, err = of$train, set = "the 25 farms it was fitted to"),
            data.frame(k = of$k, err = of$test,  set = "the other 25 farms"))
ggplot(ol, aes(k, err / 1000, colour = set)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  scale_colour_manual(values = c("blue", "red"), name = "Prediction error on") +
  expand_limits(y = 0) +
  labs(x = "Number of predictors in the model", y = "Prediction error, RMSE (t/ha)") +
  theme_classic()
of$k[which.min(of$test)]     # best on new farms: 1 predictor

# predicted vs observed on the other 25 farms, 1 vs 15 predictors
pv <- do.call(rbind, lapply(c(1, 15), function(k) {
  m <- lm(reformulate(ord[1:k], "biomass"), data = farms[tr, ])
  data.frame(model = paste(k, "predictor(s)"), obs = farms$biomass[-tr], pred = predict(m, farms[-tr, ]))
}))
ggplot(pv, aes(obs / 1000, pred / 1000)) +
  geom_abline(linetype = 2, colour = "grey50") +
  geom_point(colour = "red") +
  facet_wrap(~ model) + coord_equal() +
  labs(x = "Observed biomass, the other 25 farms (t/ha)", y = "Predicted (t/ha)") +
  theme_classic()

# Overfitting: the extra coefficients describe the noise in the first
# 25 farms, and that noise is not in the next 25.

# ---- 5) What is the model for? ---------------------------------------

# Three goals, each with its own way to choose (Tredennick et al. 2021,
# Ecology):

# | Goal | How to choose | Where | |---|---|---| | Explore: many
# variables, no strong expectations | search all subsets; the findings
# are new hypotheses | Part D | | Test ideas: you had them before the
# data | a few models written first; AIC weights | Part C | | Predict:
# new places, years, plants | held-out data, cross-validation | Part F,
# Week 6 |

# All three need a score that does not reward complexity for free. That
# is AIC (Part B).

########################################################################
# PART B. Scoring models
########################################################################

# ---- 6) Likelihood, one farm at a time -------------------------------

# Plot the raw data first, then the model on top.

# [HANDY]
# fitted values with a 95% CI, for any lm() or glm(), on the scale of y
fit_ci <- function(model, newdata, level = 0.95) {
  p <- predict(model, newdata = newdata, se.fit = TRUE)   # link scale for a glm
  q <- if (inherits(model, "glm")) qnorm(1 - (1 - level) / 2) else qt(1 - (1 - level) / 2, df.residual(model))
  inv <- family(model)$linkinv                            # identity for lm()
  fit <- as.numeric(p$fit); se <- as.numeric(p$se.fit)
  data.frame(newdata, fit = inv(fit), lwr = inv(fit - q * se), upr = inv(fit + q * se))
}
g_gdd <- data.frame(gdd = seq(min(farms$gdd), max(farms$gdd), length.out = 100))

p_gdd <- ggplot(farms, aes(gdd, biomass / 1000)) +
  geom_point(colour = "grey50") +
  labs(x = "Growing degree days (Oct to Mar)", y = "Biomass (t/ha)") +
  theme_classic()
p_gdd

# [CORE]
m_temp <- lm(biomass ~ gdd, data = farms)
logLik(m_temp)     # -420.47, df = 3: intercept, slope, and sigma

ci_temp <- fit_ci(m_temp, g_gdd)
p_gdd +
  geom_ribbon(data = ci_temp, aes(gdd, ymin = lwr / 1000, ymax = upr / 1000), inherit.aes = FALSE, fill = "blue", alpha = 0.2) +
  geom_line(data = ci_temp, aes(gdd, fit / 1000), colour = "blue", linewidth = 1)

# The model puts a normal bell around the line at each farm's degree
# days, with width sigma. A farm's likelihood is the height of its bell
# where the farm actually landed: high if the model expected it, low if
# not. The likelihood of the data is the product of the 50 heights; we
# add their logs instead. It is the probability of the data under the
# model, not the probability that the model is true. logLik() uses the
# maximum-likelihood sigma, which divides by n, not n − 2.

# [HANDY]
sigma <- sqrt(mean(resid(m_temp)^2))                  # width of the bell

# one farm: the height of the normal curve at its biomass
y <- farms$biomass[4]; y_hat <- fitted(m_temp)[[4]]
1 / (sigma * sqrt(2 * pi)) * exp(-(y - y_hat)^2 / (2 * sigma^2))
dnorm(y, mean = y_hat, sd = sigma)                    # the same, built in

# every farm at once
scores <- dnorm(farms$biomass, mean = fitted(m_temp),
                sd = sigma, log = TRUE)               # log height, each farm
round(range(scores), 2)                               # worst and best farm
sum(scores)                                           # = logLik(m_temp)

# A worse model expects the data less. The null model gives every farm
# the same prediction, so its bell is wide and low:

# [HANDY]
m_null <- lm(biomass ~ 1, data = farms)
logLik(m_null)     # -445.07, lower than -420.47

# lm() and glm() choose the coefficients that make the observed data
# most probable: maximum likelihood. Try every slope, and log L peaks
# at the slope lm() reports:

# [DEMO: log L for every slope (the peak is the lm() slope)]
slopes <- seq(0, 14, by = 0.05)
ll <- sapply(slopes, function(b) {
  a <- mean(farms$biomass) - b * mean(farms$gdd)      # best intercept for this slope
  r <- farms$biomass - a - b * farms$gdd
  sum(dnorm(r, 0, sqrt(mean(r^2)), log = TRUE))
})
slopes[which.max(ll)]     # 7.00, the lm() slope (6.99) to the grid step
plot(slopes, ll, type = "l", xlab = "Slope tried", ylab = "log L")
abline(v = coef(m_temp)["gdd"], col = "blue", lty = 2)

# For yes/no data there are no bells: each record's likelihood is the
# probability the model gave to what actually happened, p if it
# happened and 1 − p if not.

# [DEMO: likelihood of a logistic regression by hand (Wednesday's pistachios)]
pis_l <- read.csv("https://greymonroe.github.io/PLS_206/data/pistacia_niche.csv")
pis_l$y <- as.integer(pis_l$species_short == "terebinthus")
m_bin <- glm(y ~ bio1_annual_mean_temp, family = binomial, data = pis_l)
p_ter <- fitted(m_bin)                                # P(terebinthus) for each record
sum(log(ifelse(pis_l$y == 1, p_ter, 1 - p_ter)))              # = logLik(m_bin)
logLik(m_bin)

# [DEMO: the bells on the slides]
sig <- sigma
xs <- c(750, 1000, 1250)
dens <- do.call(rbind, lapply(xs, function(x0) {
  mu <- predict(m_temp, data.frame(gdd = x0))
  y <- seq(mu - 3 * sig, mu + 3 * sig, length.out = 200)
  data.frame(x0 = x0, y = y, x = x0 + 9e4 * dnorm(y, mu, sig))
}))
ggplot(farms, aes(gdd, biomass / 1000)) +
  geom_point(colour = "grey50") +
  geom_abline(intercept = coef(m_temp)[1] / 1000, slope = coef(m_temp)[2] / 1000, colour = "blue", linewidth = 1) +
  geom_path(data = dens, aes(x, y / 1000, group = x0), colour = "orange", linewidth = 1) +
  geom_vline(xintercept = xs, colour = "orange", linetype = 3) +
  labs(x = "Growing degree days (Oct to Mar)", y = "Biomass (t/ha)") +
  theme_classic()

# ---- 7) AIC = misfit + a price for each parameter --------------------

# AIC = −2 log L + 2k. Here 840.93 + 2 × 3 = 846.93. AIC() gives the
# same, as long as k counts sigma.

# [CORE]
AIC(m_temp)
-2 * as.numeric(logLik(m_temp)) + 2 * 3     # k = 3: intercept, slope, sigma
attr(logLik(m_temp), "df")                  # the k that R uses

# ---- 8) The null model -----------------------------------------------

# No predictors: k = 2 (the mean and sigma). AIC = −2 × −445.07 + 4 =
# 894.1, about 47 worse than the gdd model.

# [CORE]
m_null <- lm(biomass ~ 1, data = farms)
logLik(m_null)
AIC(m_null)
AIC(m_null) - AIC(m_temp)

ci_null <- fit_ci(m_null, g_gdd)
p_gdd +
  geom_ribbon(data = ci_null, aes(gdd, ymin = lwr / 1000, ymax = upr / 1000), inherit.aes = FALSE, fill = "grey40", alpha = 0.2) +
  geom_line(data = ci_null, aes(gdd, fit / 1000), colour = "grey30", linewidth = 1) +
  geom_line(data = ci_temp, aes(gdd, fit / 1000), colour = "blue", linewidth = 1)

# ---- 9) Misfit falls, the penalty rises, AIC finds the balance -------

# All 50 farms, adding predictors best first. The misfit drops sharply
# when gdd comes in, then barely moves; the penalty rises by 2 with
# every predictor. AIC is lowest at 2 predictors.

# [DEMO: misfit, penalty and AIC as predictors are added]
ord50 <- names(sort(-abs(cor(farms)[P, "biomass"])))     # best first, all 50 farms
fp <- do.call(rbind, lapply(0:15, function(k) {
  m <- lm(if (k == 0) biomass ~ 1 else reformulate(ord50[1:k], "biomass"), data = farms)
  h <- hatvalues(m)
  data.frame(k = k, misfit = -2 * as.numeric(logLik(m)), penalty = 2 * attr(logLik(m), "df"),
             aic = AIC(m), bic = BIC(m),
             loo = sqrt(mean((resid(m) / (1 - h))^2)))   # leave-one-out error, used in Part F
}))
fp$k[which.min(fp$aic)]    # 2 predictors
fl <- rbind(data.frame(k = fp$k, v = fp$misfit - min(fp$misfit), part = "misfit: -2 log L"),
            data.frame(k = fp$k, v = fp$penalty - min(fp$penalty), part = "penalty: 2k"),
            data.frame(k = fp$k, v = fp$aic - min(fp$aic), part = "AIC = misfit + penalty"))
fl$part <- factor(fl$part, levels = unique(fl$part))
ggplot(fl, aes(k, v, colour = part)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  scale_colour_manual(values = c("blue", "red", "black"), name = NULL) +
  labs(x = "Number of predictors in the model", y = "Change from its lowest value") +
  theme_classic()

# ---- 10) ΔAIC and Akaike weights -------------------------------------

# Only differences mean anything: ΔAIC = AIC − the lowest AIC in the
# set. Rules of thumb (Burnham and Anderson): within about 2, about
# equally supported; 4 to 7, clearly less; over 10, essentially none.
# Each model's weight is exp(−Δ/2), divided by the total, so the
# weights add to 1. Here are three of the models from Part C (m_tm adds
# rain, m_all adds rain and soil):

# [DEMO: delta and weights by hand]
m_tm  <- lm(biomass ~ gdd + rain_total, data = farms)
m_all <- lm(biomass ~ gdd + rain_total + organic_matter + soil_ph, data = farms)
a <- c(m_temp = AICc(m_temp), m_tm = AICc(m_tm), m_all = AICc(m_all))
delta <- a - min(a)
w <- exp(-delta / 2) / sum(exp(-delta / 2))
round(cbind(AICc = a, delta, weight = w), 2)    # deltas 0, 2.11, 5.89; weights 0.71, 0.25, 0.04

# ---- 11) AICc: a correction for small samples ------------------------

# AICc = AIC + 2k(k + 1) / (n − k − 1). Large when k is close to n,
# fading as n grows. Use it when n/k is below about 40; with 50 farms
# it adds about half a point.

# [CORE]
AICc(m_temp)     # MuMIn

# [DEMO: AICc by hand]
n <- nobs(m_temp); k <- 3
AIC(m_temp) + 2 * k * (k + 1) / (n - k - 1)

# ---- 12) BIC: a stiffer price ----------------------------------------

# BIC = −2 log L + k log(n). With 50 farms, log(50) = 3.91 per
# parameter instead of 2; with the 6,299 Pistacia records in Part D,
# 8.7. BIC picks smaller models than AIC, and the difference grows with
# n. AIC aims at the model that predicts new data best; BIC at the true
# model, if it is in the set.

# [CORE]
BIC(m_temp)
-2 * as.numeric(logLik(m_temp)) + 3 * log(50)
log(c(50, 6299))   # BIC's price per parameter; AIC's is 2

# [DEMO: the penalty per parameter as n grows]
pen <- data.frame(n = exp(seq(log(5), log(20000), length.out = 300)))
ggplot(pen, aes(n, log(n))) +
  geom_line(colour = "purple", linewidth = 1) +
  geom_hline(yintercept = 2, colour = "blue", linewidth = 1) +
  annotate("point", x = c(50, 6299), y = log(c(50, 6299)), size = 3) +
  annotate("text", x = 3000, y = 1.5, label = "AIC: 2 per parameter", colour = "blue") +
  annotate("text", x = 200, y = 7.5, label = "BIC: log(n) per parameter", colour = "purple") +
  scale_x_log10() +
  labs(x = "Sample size, n (log scale)", y = "Penalty per parameter") +
  theme_classic()

########################################################################
# PART C. Testing ideas
########################################################################

# ---- 13) Ideas first, then models: model.sel() -----------------------

# Three ideas about what limits cover-crop growth over winter, each
# written as a model before looking: temperature (gdd), moisture
# (rain_total), soil (organic_matter + soil_ph). Add the combinations
# that make sense and the null. model.sel() fits nothing itself: it
# ranks the models you give it by AICc, with delta and weights.

# [CORE]
m_null  <- lm(biomass ~ 1, data = farms)
m_temp  <- lm(biomass ~ gdd, data = farms)
m_moist <- lm(biomass ~ rain_total, data = farms)
m_soil  <- lm(biomass ~ organic_matter + soil_ph, data = farms)
m_tm    <- lm(biomass ~ gdd + rain_total, data = farms)
m_all   <- lm(biomass ~ gdd + rain_total + organic_matter + soil_ph, data = farms)
ms <- model.sel(m_null, m_temp, m_moist, m_soil, m_tm, m_all)
round(as.data.frame(ms)[, c("df", "logLik", "AICc", "delta", "weight")], 2)

# Temperature alone is best (weight 0.71). Temperature + moisture is
# 2.11 behind with almost the same log likelihood: the gap is the
# penalty for one more parameter. Moisture alone and soil alone are
# worse than the null.

# [HANDY]
ms                                 # the full table: coefficients of every model
model.sel(m_null, m_temp, m_moist, m_soil, m_tm, m_all, rank = BIC)   # ranked by BIC instead

# ---- 14) The evidence, as weights ------------------------------------

# [HANDY]
wd <- data.frame(model = rownames(ms), weight = ms$weight)
wd$model <- factor(wd$model, levels = rev(wd$model))
ggplot(wd, aes(weight, model)) +
  geom_col(fill = "blue", width = 0.6) +
  geom_text(aes(label = sprintf("%.2f", weight)), hjust = -0.2) +
  scale_x_continuous(limits = c(0, 1)) +
  labs(x = "Akaike weight", y = NULL) +
  theme_classic()

# ---- 15) Nested models: anova() and AIC agree ------------------------

# m_tm is m_temp plus one term, so the F test from Week 2 works too: p
# = 0.62 and ΔAICc = 2.11, the same conclusion.

# [CORE]
anova(m_temp, m_tm)      # does adding moisture help?

# ---- 16) Not nested? Only AIC can compare them -----------------------

# Moisture vs soil: neither model is the other plus extra terms.
# anova() runs, but gives a negative sum of squares and no test. AIC
# can compare any models fitted to the same rows with the same
# response.

# [CORE]
anova(m_moist, m_soil)   # moisture vs soil: not nested
AICc(m_moist) - AICc(m_soil)

# ---- 17) Some terms stay, no matter what -----------------------------

# Not up for selection: design terms (blocks, plots: the analysis must
# match how the data were collected), known controls (sampling effort
# in Part D, kept with fixed =), and the main effects under any
# interaction or the linear term under a quadratic (dredge() respects
# this by default).

# ---- 18) Same response: log(y) and y cannot be compared --------------

# [CORE]
m_raw <- lm(biomass ~ gdd, data = farms)
m_log <- lm(log(biomass) ~ gdd, data = farms)
AIC(m_raw, m_log)        # 846.9 vs 123.1: different units, not a better model

# The likelihood of log biomass is in different units from the
# likelihood of biomass. Pick the scale first, from the residuals (Week
# 2), then compare models on that scale.

# ---- 19) Same rows: missing values change n --------------------------

# Suppose three soil tests were lost. lm() quietly drops those farms
# from any model that uses soil_ph, so that model is scored on 47 farms
# and the other on 50.

# [CORE]
farms_na <- farms
farms_na$soil_ph[c(3, 17, 29)] <- NA    # suppose 3 soil tests were lost
m_a <- lm(biomass ~ gdd, data = farms_na)
m_b <- lm(biomass ~ gdd + soil_ph, data = farms_na)
c(nobs(m_a), nobs(m_b))
AIC(m_a, m_b)            # m_b looks 48 better only because it has 3 fewer farms

# The fix: drop the incomplete rows once, before fitting any model.

# [CORE]
farms_cc <- na.omit(farms_na)
m_a <- lm(biomass ~ gdd, data = farms_cc)
m_b <- lm(biomass ~ gdd + soil_ph, data = farms_cc)
c(nobs(m_a), nobs(m_b))
AIC(m_a, m_b)

# ---- 20) Report the table, not just the winner -----------------------

# Report the candidate set and the table (k, AICc, ΔAICc, weight), then
# the best model the usual way: coefficients with confidence intervals
# and a plot of the fit with the raw data (section 6).

# [CORE]
round(as.data.frame(ms)[, c("df", "AICc", "delta", "weight")], 2)
summary(m_temp)
confint(m_temp)

########################################################################
# PART D. Exploring
########################################################################

# ---- 21) The data: two Mediterranean pistachios ----------------------

# Pistacia lentiscus (mastic tree, evergreen) and P. terebinthus
# (turpentine tree, deciduous), wild relatives of the pistachio around
# the Mediterranean. 6,299 GBIF records with the 19 WorldClim
# bioclimatic variables at each spot. The response is y = 1 for
# terebinthus, 0 for lentiscus: a logistic regression (Week 3).

# [CORE]
pis <- read.csv("https://greymonroe.github.io/PLS_206/data/pistacia_niche.csv")
names(pis) <- sub("^(bio[0-9]+)_.*", "\\1", names(pis))   # bio1 ... bio19
pis$y <- as.integer(pis$species_short == "terebinthus")
table(pis$species_short)
length(unique(pis$country))

# [DEMO: the map (needs the maps package)]
# install.packages("maps")
if (requireNamespace("maps", quietly = TRUE)) {
  world <- map_data("world")
  print(ggplot() +
    geom_polygon(data = world, aes(long, lat, group = group), fill = "grey93", colour = "grey75", linewidth = 0.2) +
    geom_point(data = pis, aes(lon, lat, colour = species_short), size = 0.5, alpha = 0.5) +
    scale_colour_manual(values = c(lentiscus = "darkgreen", terebinthus = "sienna"), name = NULL) +
    coord_quickmap(xlim = c(-12, 44), ylim = c(28, 47)) +
    theme_void())
}

# ---- 22) 19 climate variables, all tangled ---------------------------

# The species differ, but overlap a lot, and the 19 variables are built
# from the same monthly temperatures and rainfall.

# [HANDY]
bios <- paste0("bio", 1:19)
v4 <- c(bio1 = "Annual mean temp (°C)", bio6 = "Coldest-month min (°C)", bio4 = "Temp seasonality", bio12 = "Annual rain (mm)")
dl <- do.call(rbind, lapply(names(v4), function(v) data.frame(var = v4[[v]], x = pis[[v]], species = pis$species_short)))
dl$var <- factor(dl$var, levels = v4)
ggplot(dl, aes(x, fill = species, colour = species)) +
  geom_density(alpha = 0.3) +
  facet_wrap(~ var, scales = "free", nrow = 1) +
  scale_fill_manual(values = c("darkgreen", "sienna")) + scale_colour_manual(values = c("darkgreen", "sienna")) +
  labs(x = NULL) +
  theme_classic() + theme(legend.position = "bottom")
round(cor(pis[bios])[1:6, 1:6], 2)

# A first model: species ~ temperature. Before choosing among 19
# variables, one logistic regression (Week 3). Each record's likelihood
# is the probability the model gave to what happened: the curve's
# height for a terebinthus, 1 minus it for a lentiscus.

# [HANDY]
m_1 <- glm(y ~ bio1, family = binomial, data = pis)
p_ter <- fitted(m_1)                                   # P(terebinthus), each record
sum(log(ifelse(pis$y == 1, p_ter, 1 - p_ter)))         # likelihood by hand
logLik(m_1)                                            # the same

# ---- 23) Throw everything in? bio7 is bio5 − bio6 --------------------

# All 19 in one logistic regression: three coefficients near ±30,000,
# for variables measured in degrees. bio7 (annual temperature range) is
# bio5 (warmest-month maximum) − bio6 (coldest-month minimum), so no
# data can separate the three.

# [CORE]
m_all <- glm(y ~ ., family = binomial,
             data = pis[, c("y", paste0("bio", 1:19))])
round(coef(m_all)[c("bio5", "bio6", "bio7")])
range(pis$bio7 - (pis$bio5 - pis$bio6))   # bio7 IS bio5 - bio6

# Before going on: the sampling-effort control (section 26) on a log
# scale, the records that have it, and the climate variables in SD
# units (scale()) so their coefficients compare. Keep the means and SDs
# to put the results back in raw units (section 32).

# [CORE]
pis$effort <- log10(pis$plant_obs_effort)
pis <- subset(pis, !is.na(effort))          # the same rows for every model
nrow(pis)
sc <- scale(pis[bios])
bio_mean <- attr(sc, "scaled:center"); bio_sd <- attr(sc, "scaled:scale")
pis[bios] <- sc

# ---- 24) No pair is correlated, yet one is redundant -----------------

# A popular shortcut: keep variables, best single predictor first, so
# that no pair correlates above 0.7. It keeps 9 variables, all pairs
# below 0.66, but VIFs (Week 2) up to 79: a variable can be almost
# perfectly predicted by a combination of the others.

# [DEMO: the pairwise filter]
aics <- sapply(setdiff(bios, "bio7"), function(v) AIC(glm(reformulate(v, "y"), family = binomial, data = pis)))
cc <- cor(pis[bios])
pw <- c()
for (v in names(sort(aics))) if (all(abs(cc[v, pw]) < 0.7)) pw <- c(pw, v)
pw

# [CORE]
max(abs(cor(pis[pw])[upper.tri(diag(length(pw)))]))
round(vif(glm(reformulate(pw, "y"), family = binomial, data = pis)), 1)

# What a VIF is: regress one predictor on all the others; VIF = 1 / (1
# − R²). bio6 is predicted almost perfectly by the other eight
# together, though no single pair is above 0.7.

# [HANDY]
r2 <- summary(lm(reformulate(setdiff(pw, "bio6"), "bio6"), data = pis))$r.squared
r2              # 0.989
1 / (1 - r2)    # its VIF (car::vif() on the glm weights by the fit, so it gives a bit less)

# ---- 25) Drop the worst VIF, refit, repeat ---------------------------

# The Week 2 loop: fit all the variables, drop the one with the highest
# VIF, refit, and repeat until every VIF is under 5. 11 rounds, from 18
# variables to 7. Which variable of a correlated group survives is
# partly luck: say which you kept and why.

# [CORE]
keep <- setdiff(bios, "bio7")
repeat {
  v <- vif(glm(reformulate(keep, "y"), family = binomial, data = pis))
  if (max(v) < 5) break
  keep <- setdiff(keep, names(which.max(v)))   # drop the worst, refit
}
keep
round(v, 1)

# [DEMO: the largest VIF at each round]
k_vars <- setdiff(bios, "bio7"); steps <- NULL
repeat {
  v <- vif(glm(reformulate(k_vars, "y"), family = binomial, data = pis))
  steps <- rbind(steps, data.frame(n = length(k_vars), maxvif = max(v), worst = names(which.max(v))))
  if (max(v) < 5) break
  k_vars <- setdiff(k_vars, names(which.max(v)))
}
steps
ggplot(steps, aes(n, maxvif)) +
  geom_hline(yintercept = 5, colour = "red", linetype = 2) +
  geom_line() + geom_point() +
  geom_text(aes(label = worst), vjust = -0.8, size = 3) +
  scale_x_reverse() + scale_y_log10() +
  labs(x = "Climate variables still in the model", y = "Largest VIF (log scale)") +
  theme_classic()

# ---- 26) Sampling effort: a control we keep --------------------------

# Some places are surveyed far more than others. effort is the number
# of plant records of any species in the same pixel, log10. We know
# sampling effort matters for records like these, so it goes in every
# model, chosen before looking.

# [HANDY]
p_eff <- ggplot(pis, aes(effort, y)) +
  geom_jitter(height = 0.03, width = 0, alpha = 0.05, colour = "grey40") +
  labs(x = "Plant records in the same pixel, all species (log10)", y = "P(P. terebinthus)") +
  theme_classic()
p_eff

# [CORE]
m_eff <- glm(y ~ effort, family = binomial, data = pis)
round(summary(m_eff)$coefficients, 3)

ci_eff <- fit_ci(m_eff, data.frame(effort = seq(min(pis$effort), max(pis$effort), length.out = 200)))
p_eff +
  geom_ribbon(data = ci_eff, aes(effort, ymin = lwr, ymax = upr), inherit.aes = FALSE, fill = "blue", alpha = 0.25) +
  geom_line(data = ci_eff, aes(effort, fit), colour = "blue", linewidth = 1)

# ---- 27) dredge(): fit every subset ----------------------------------

# dredge() fits every subset of the predictors in a full model and
# ranks them by AICc. It refuses to run unless missing values cause an
# error (the same-rows rule, enforced), and fixed = ~ effort keeps
# effort in every model: 2^7 = 128 models.

# [CORE]
options(na.action = "na.fail")     # dredge() insists
full <- glm(y ~ bio1 + bio3 + bio4 + bio8 + bio9 + bio18 + bio19 + effort,
            family = binomial, data = pis)
dd <- dredge(full, fixed = ~ effort)   # effort in every model
nrow(dd)

# ---- 28) The top 10 of 128 -------------------------------------------

# The top model has six of the seven (not bio9), with only 37% of the
# weight, and 3 models are within 2 AICc of it. No single model is the
# answer. (Printing a few rows of a dredge() table rescales their
# weights to add to 1; recalc.weights = FALSE keeps the weights from
# all 128.)

# [CORE]
dd[1:10, recalc.weights = FALSE]                  # keep the weights from all 128
subset(dd, delta < 2, recalc.weights = FALSE)

# [DEMO: the tile picture, and all 128 models]
x <- as.data.frame(dd)[1:10, ]
tiles <- do.call(rbind, lapply(1:10, function(i) data.frame(rank = i, var = keep, inm = !is.na(unlist(x[i, keep])))))
tiles$var <- factor(tiles$var, levels = keep)
ggplot(tiles, aes(var, rank, fill = inm)) +
  geom_tile(colour = "white", linewidth = 1) +
  scale_fill_manual(values = c(`TRUE` = "blue", `FALSE` = "grey92"), guide = "none") +
  scale_y_reverse(breaks = 1:10, labels = sprintf("%d   (delta %.2f, w %.2f)", 1:10, x$delta, x$weight)) +
  scale_x_discrete(position = "top") +
  labs(x = NULL, y = NULL) +
  theme_classic()

ggplot(data.frame(rank = seq_len(nrow(dd)), delta = dd$delta), aes(rank, delta)) +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = 2, fill = "green", alpha = 0.2) +
  geom_point() +
  scale_y_sqrt() +
  labs(x = "Model rank (all 128 models)", y = "Delta AICc (square-root scale)") +
  theme_classic()

# ---- 29) AICc vs BIC with 6,299 records ------------------------------

# Ranked by BIC, the top model drops bio3 and bio18: 6 parameters
# instead of 8. BIC charges log(6299) = 8.7 per parameter, more than
# four times AIC's 2. With large samples, AIC keeps small effects
# because it has the data to estimate them; BIC drops them. Say which
# one you used.

# [CORE]
ddb <- dredge(full, fixed = ~ effort, rank = "BIC")
ddb[1:5, recalc.weights = FALSE]
dd[1:5, recalc.weights = FALSE]

# ---- 30) Importance: summed weights ----------------------------------

# For each variable, add up the weights of the models that contain it.
# With seven variables in or out, every variable is in exactly half the
# models, so even a useless one can collect some weight. Importance is
# not an effect size and not a p-value.

# [CORE]
sw(dd)        # importance: summed Akaike weights

# ---- 31) Model averaging: estimates across the top models ------------

# Average the coefficients over the well-supported models, weighted by
# their Akaike weights. The full average counts a zero for a model
# without the variable; the conditional average uses only the models
# that contain it. Variables in every top model (bio1, bio4) barely
# change; variables in only some (bio3, bio9) are pulled toward zero.

# [CORE]
avg <- model.avg(get.models(dd, subset = delta < 2))
summary(avg)
round(coef(avg, full = TRUE), 2)

# [DEMO: averaged vs top-model coefficients]
ct <- coefTable(avg, full = TRUE)
top <- get.models(dd, 1)[[1]]
tt <- summary(top)$coefficients
cl <- rbind(data.frame(var = rownames(ct), est = ct[, 1], se = ct[, 2], what = "averaged (delta < 2)"),
            data.frame(var = rownames(tt), est = tt[, 1], se = tt[, 2], what = "top model only"))
cl <- subset(cl, var %in% keep)
cl$var <- factor(cl$var, levels = rev(keep))
ggplot(cl, aes(est, var, colour = what)) +
  geom_vline(xintercept = 0, colour = "grey60") +
  geom_pointrange(aes(xmin = est - 1.96 * se, xmax = est + 1.96 * se), position = position_dodge(width = 0.5)) +
  scale_colour_manual(values = c("blue", "grey45"), name = NULL) +
  labs(x = "Coefficient (log odds of P. terebinthus per SD)", y = NULL) +
  theme_classic() + theme(legend.position = "bottom")

# ---- 32) What we found: hypotheses, not results ----------------------

# The pattern from the averaged model, one variable at a time: predict
# over a grid of that variable with the other climate variables at
# their mean (0 in SD units) and effort at its mean, then put the
# x-axis back in raw units. The raw records go along the top and
# bottom.

# [CORE]
partial_avg <- function(v) {
  nd <- data.frame(matrix(0, 150, length(keep), dimnames = list(NULL, keep)))
  nd$effort <- mean(pis$effort)
  nd[[v]] <- seq(quantile(pis[[v]], 0.01), quantile(pis[[v]], 0.99), length.out = 150)
  pr <- predict(avg, nd, se.fit = TRUE, full = TRUE)          # link (log-odds) scale
  data.frame(var = v, x = nd[[v]] * bio_sd[[v]] + bio_mean[[v]],
             fit = plogis(pr$fit), lwr = plogis(pr$fit - 1.96 * pr$se.fit), upr = plogis(pr$fit + 1.96 * pr$se.fit))
}
pe <- rbind(partial_avg("bio4"), partial_avg("bio1"))
raw_pts <- rbind(data.frame(var = "bio4", x = pis$bio4 * bio_sd[["bio4"]] + bio_mean[["bio4"]], y = pis$y),
                 data.frame(var = "bio1", x = pis$bio1 * bio_sd[["bio1"]] + bio_mean[["bio1"]], y = pis$y))
labs_v <- c(bio1 = "Annual mean temperature (bio1, °C)", bio4 = "Temperature seasonality (bio4, SD x 100)")

p_raw <- ggplot(raw_pts, aes(x, y)) +
  geom_jitter(height = 0.03, width = 0, alpha = 0.04, size = 0.6, colour = "grey40") +
  facet_wrap(~ var, scales = "free_x", labeller = as_labeller(labs_v)) +
  labs(x = NULL, y = "P(P. terebinthus)") +
  theme_classic()
p_raw
p_raw +
  geom_ribbon(data = pe, aes(x, ymin = lwr, ymax = upr), inherit.aes = FALSE, fill = "blue", alpha = 0.25) +
  geom_line(data = pe, aes(x, fit), colour = "blue", linewidth = 1)

# Terebinthus is more likely in cooler places with stronger seasons,
# which fits it being the deciduous one of the pair. But we searched
# 128 models for this: it is a hypothesis to test with new data, not a
# result.

########################################################################
# PART E. Traps
########################################################################

# ---- 33) Stepwise selection: drop one at a time ----------------------

# stepAIC() (MASS) starts from the full model and drops or adds one
# variable at a time while AIC improves. It visits a tiny fraction of
# the 32,768 models. On the cover-crop farms it keeps gdd and
# tmean_winter; the truth is gdd alone.

# [HANDY]
full_f <- lm(biomass ~ ., data = farms)
step_f <- stepAIC(full_f, trace = 0)     # MASS
step_f$anova                             # the path: one predictor dropped per step
attr(terms(step_f), "term.labels")       # what it kept

# forward: start with none, add the best one each step
step_fw <- stepAIC(lm(biomass ~ 1, data = farms), direction = "forward", trace = 0,
                   scope = list(lower = ~ 1, upper = reformulate(setdiff(names(farms), "biomass"))))
step_fw$anova                         # + gdd, + tmean_winter, stop

# ---- 34) Resample the farms, rerun stepAIC ---------------------------

# Resample the 50 farms with replacement 200 times (as if collecting
# the data again) and run stepAIC() each time. gdd is kept 99% of the
# time, but every one of the 14 variables that do nothing is also kept
# in 35% to 57.5% of the runs. The number kept ranges from 2 to 13
# (median 8), and not one run finds gdd alone. (About a minute.)

# [DEMO: stepAIC on 200 resampled data sets]
set.seed(2)
sel <- replicate(200, {
  s <- farms[sample(50, 50, replace = TRUE), ]
  m <- stepAIC(lm(biomass ~ ., data = s), trace = 0)
  P %in% attr(terms(m), "term.labels")
})
rownames(sel) <- P
freq <- sort(rowMeans(sel), decreasing = TRUE)
round(freq, 3)                          # share of runs that kept each variable
nk <- colSums(sel)                      # predictors kept in each run
summary(nk)
sum(nk == 1 & sel["gdd", ])             # runs that found the truth: gdd alone

fr <- data.frame(var = factor(names(freq), levels = rev(names(freq))), f = freq)
ggplot(fr, aes(f, var, fill = var == "gdd")) +
  geom_col(width = 0.7) +
  scale_fill_manual(values = c(`TRUE` = "blue", `FALSE` = "grey60"), guide = "none") +
  scale_x_continuous(limits = c(0, 1), labels = scales::percent) +
  labs(x = "Share of 200 resampled data sets where stepAIC() kept it", y = NULL) +
  theme_classic()
ggplot(data.frame(k = nk), aes(k)) +
  geom_bar(fill = "grey60") +
  geom_vline(xintercept = 1, colour = "blue", linewidth = 1) +
  scale_x_continuous(breaks = 0:15, limits = c(0, 15.5)) +
  labs(x = "Predictors stepAIC() kept (the truth: 1)", y = "Resampled data sets") +
  theme_classic()

# ---- 35) Pure noise, 'significant' results ---------------------------

# 100 rows, a response and 20 predictors, all random numbers. Stepwise
# keeps 5 of the 20, two of them with p = 0.006 and 0.001. The p-values
# are wrong because the same data were used twice: once to choose the
# variables, and again to test them.

# [DEMO: stepAIC on pure noise]
set.seed(41)
noise <- data.frame(y = rnorm(100), matrix(rnorm(100 * 20), 100))
step_n <- stepAIC(lm(y ~ ., data = noise), trace = 0)
round(summary(step_n)$coefficients, 3)

# ---- 36) Fix: choose on one half, test on the other ------------------

# 1,000 pure-noise data sets of 100 rows. Choose with stepAIC() on the
# first 50 rows and test on the same 50: 86% end with at least one
# predictor at p < 0.05. Test the chosen variables on the other 50 rows
# instead: 25%, about what chance gives for testing five or six
# variables at once. (About a minute.)

# [DEMO: 1,000 pure-noise data sets]
set.seed(1)
fz <- t(replicate(1000, {
  nz <- data.frame(y = rnorm(100), matrix(rnorm(100 * 20), 100))
  a <- nz[1:50, ]; b <- nz[51:100, ]                 # choose on a, test on b
  m <- stepAIC(lm(y ~ ., data = a), trace = 0)
  tl <- attr(terms(m), "term.labels")
  same <- if (length(tl)) any(summary(m)$coefficients[-1, 4] < 0.05) else FALSE
  conf <- if (length(tl)) any(summary(lm(reformulate(tl, "y"), data = b))$coefficients[-1, 4] < 0.05) else FALSE
  c(same = same, conf = conf, k = length(tl))
}))
100 * colMeans(fz[, c("same", "conf")])   # % with a 'significant' predictor
mean(fz[, "k"])                           # predictors kept, on average

nb <- data.frame(what = c("chosen and tested\non the same 50 rows", "chosen on 50 rows,\ntested on the other 50"),
                 pct = colMeans(fz[, c("same", "conf")]))
nb$what <- factor(nb$what, levels = nb$what)
ggplot(nb, aes(what, pct)) +
  geom_col(fill = c("red", "darkgreen"), width = 0.5) +
  geom_text(aes(label = scales::percent(pct, accuracy = 1)), vjust = -0.4) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
  labs(x = NULL, y = "Data sets with a 'significant' predictor") +
  theme_classic()

# ---- 37) Selection finds predictors, not causes ----------------------

# Latitude alone beats the null by 6.9 AICc (97% of the weight). In the
# simulation latitude does nothing directly: farms further north have
# fewer degree days, and degree days drive biomass. AIC rewards
# whatever predicts, not what causes.

# [CORE]
m_lat <- lm(biomass ~ latitude, data = farms)
round(as.data.frame(model.sel(m_null, m_lat))[, c("df", "AICc", "delta", "weight")], 3)

g_lat <- data.frame(latitude = seq(min(farms$latitude), max(farms$latitude), length.out = 100))
ci_lat <- fit_ci(m_lat, g_lat)
p_lat <- ggplot(farms, aes(latitude, biomass / 1000)) +
  geom_point(colour = "grey50") +
  labs(x = "Latitude (°N)", y = "Biomass (t/ha)") +
  theme_classic()
p_lat
p_lat +
  geom_ribbon(data = ci_lat, aes(latitude, ymin = lwr / 1000, ymax = upr / 1000), inherit.aes = FALSE, fill = "blue", alpha = 0.2) +
  geom_line(data = ci_lat, aes(latitude, fit / 1000), colour = "blue", linewidth = 1)

# [HANDY]
cor(farms$latitude, farms$gdd)     # the path: latitude -> degree days -> biomass

########################################################################
# PART F. Predicting
########################################################################

# ---- 38) AIC approximates leave-one-out prediction -------------------

# AIC was designed to estimate how well a model will predict new data;
# for large samples it behaves like leave-one-out cross-validation (fit
# to all but one farm, predict that farm, repeat for every farm). For a
# linear model the leave-one-out residual is resid / (1 - hatvalue), no
# refitting needed (computed in section 9). On the cover-crop sequence
# both curves bottom out at 2 predictors.

# [DEMO: AIC vs leave-one-out error]
fp$k[which.min(fp$loo)]
al <- rbind(data.frame(k = fp$k, v = (fp$aic - min(fp$aic)) / diff(range(fp$aic)), what = "AIC"),
            data.frame(k = fp$k, v = (fp$loo - min(fp$loo)) / diff(range(fp$loo)), what = "leave-one-out prediction error"))
ggplot(al, aes(k, v, colour = what)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  scale_colour_manual(values = c("black", "red"), name = NULL) +
  labs(x = "Number of predictors (cover-crop farms)", y = "Scaled to 0-1") +
  theme_classic()

# ---- 39) Next week: shrink instead of choose -------------------------

# Selection treats every predictor as fully in or fully out. Week 6
# keeps them all and shrinks their slopes toward zero (ridge and
# lasso), choosing how hard to shrink by cross-validation. The lasso
# shrinks some slopes to exactly zero, so it selects too.

# The week on one line each:

# | Goal | Do | |---|---| | Test ideas | write the models first;
# model.sel(); report the table | | Explore | dredge(); the results are
# hypotheses | | Predict | held-out data, cross-validation | | Always |
# the same rows and the same response; keep design terms and controls |
# | Avoid | stepwise; p-values after selection |